package ua.yandex.shad.tempseries;

import java.util.Arrays;

class IllegalTemperatureValueException extends IllegalArgumentException {
    public IllegalTemperatureValueException() {}
    public IllegalTemperatureValueException(String msg) {
        super(msg);
    }
}

class NotEnoughValuesToForecastException extends IllegalArgumentException {
    public NotEnoughValuesToForecastException() {}
    public NotEnoughValuesToForecastException(String msg) {
        super(msg);
    }
}
public class TemperatureSeriesAnalysis {

    public static final int SSA_LENGTH = 3;
    public static final double ABSOLUTE_ZERO = - 273.0;
    private static final double THREE = 3.0;
    private static final double ROOT = Math.sqrt(THREE);
    private static final double EPSILON = 1e-05;
    private static final int MIN_FORECAST_SIZE = 10;
    private final double[] temperatureSeries;

    public TemperatureSeriesAnalysis(double[] temperatureSeries) {
        this.temperatureSeries = temperatureSeries;
        for (double d: temperatureSeries) {
            if(d < ABSOLUTE_ZERO) {
                throw new IllegalTemperatureValueException("Temperature lower than absolute zero "+d);
            }
        }
    }

    public double calcAverageTemp() {

        if (temperatureSeries == null || temperatureSeries.length == 0) {
            return 0;
        }
        double sum = 0;
        for (double d : temperatureSeries) {
            sum += d;
        }
        return sum / temperatureSeries.length;
    }

    public double calcTempDeviation() {
        if (temperatureSeries == null || temperatureSeries.length == 0) {
            return 0;
        }
        double mean = calcAverageTemp();
        double sum = 0;
        for (double d: temperatureSeries){
            sum += (d - mean) * (d - mean);
        }
        return sum / temperatureSeries.length;
    }

    public double findTempClosestToZero() {
        if (temperatureSeries == null || temperatureSeries.length == 0) {
            return 0;
        }
        double close = temperatureSeries[0];
        for (double d: temperatureSeries) {
            if (Math.abs(d) < Math.abs(close)) {
                close = d;
                continue;
            }
            if (Math.abs(d) == Math.abs(close) && d > 0 && close < 0) {
                close = d;
            }
        }
        return close;
    }

    public double findTempClosestToValue(final double tempValue) {
        if (temperatureSeries == null || temperatureSeries.length == 0) {
            return 0;
        }
        double close = temperatureSeries[0];
        for (double d: temperatureSeries) {
            if (Math.abs(tempValue - d) < Math.abs(tempValue - close)) {
                close = d;
                continue;
            }
            if ((Math.abs(tempValue - d) == Math.abs(tempValue - close)) &&
                (d > close)) {
                close = d;
            }
        }
        return close;
    }

    public double[] findTempsLessThen(final double tempValue) {
        if (temperatureSeries == null || temperatureSeries.length == 0) {
            return new double[]{};
        }
        int length = 0;
        for (double d : temperatureSeries) {
            if (d  < tempValue) {
                length++;
            }
        }
        double[] result = new double[length];
        int i = 0;
        for (double d: temperatureSeries) {
            if (d  < tempValue) {
                result[i] = d;
                i++;
            }
        }
        return result;
    }

    public double[] findTempsGreaterThen(final double tempValue) {
        if (temperatureSeries == null || temperatureSeries.length == 0) {
            return new double[]{};
        }
        int length = 0;
        for (double d : temperatureSeries) {
            if (d  > tempValue) {
                length++;
            }
        }
        double[] result = new double[length];
        int i = 0;
        for (double d: temperatureSeries) {
            if (d  > tempValue) {
                result[i] = d;
                i++;
            }
        }
        return result;
    }

    /**
    *Matrix N - SSA_LENGTH x SSA_LENGTH
    */
    private double[][] getHankelMatrix() {
        int SVD_LENGTH = temperatureSeries.length - SSA_LENGTH + 1;
        double[][] result = new double[SVD_LENGTH][SSA_LENGTH];
        for (int i = 0; i < SSA_LENGTH;i++) {
            for (int j = 0; j < SVD_LENGTH; j++) {
                result[j][i] = temperatureSeries[j+i];
            }
        }
        return result;
    }

    private double[][] getSVDMatrix() {
        double[][] hankelMatrix = getHankelMatrix();
        double[][] result = new double[SSA_LENGTH][SSA_LENGTH];
        int SVD_LENGTH = temperatureSeries.length -SSA_LENGTH + 1;
        for (int i = 0; i< SSA_LENGTH; i++) {
            for (int j = 0; j < SSA_LENGTH; j++) {
                result[i][j] = 0;
                for (int k = 0; k < SVD_LENGTH; k++) {
                    result[i][j] += hankelMatrix[k][i] * hankelMatrix[k][j];
                }
            }
        }
        return result;
    }

    public static double[] computeCrossProduct(double[] arg1, double[] arg2) {
       double[] result = new double[SSA_LENGTH];
       result[0] = - arg2[1] * arg1[2] + arg1[1] * arg2[2];
       result[1] = arg2[0] * arg1[2] - arg1[0] * arg2[2];
       result[2] = - arg2[0] * arg1[1] + arg1[0] * arg2[1];
       return result;
    }

    private static double getVectorNorm(double[] arg) {
        double result = 0;
        for (double d : arg) {
            result += d * d;
        }
        result = Math.sqrt(result);
        return result;
    }
    private static double[]  getFirstComplement(double[] arg1, double[] arg2) {
        double[] result = computeCrossProduct(arg2, arg1);
        double norm = getVectorNorm(result);
        for (double d: result) {
            if (norm != 0) {
                d = d / norm;
            }
        }
        return result;
    }
    private static void getSecondComplement(double[] arg1, double[] arg2, double[] arg3) {
        double[] arg1Norm = new double[SSA_LENGTH];
        for (int i = 0; i < SSA_LENGTH; i++) {
            arg1Norm[i] = arg1[i] / getVectorNorm(arg1);
        }
        if (Math.abs(arg1Norm[0]) >= Math.abs(arg1Norm[1])) {
            double inverseLength = 1 / Math.sqrt(arg1Norm[0] * arg1Norm[0] + arg1Norm[2] * arg1Norm[2]);
            arg2[0] = - arg1Norm[2] * inverseLength;
            arg2[1] = 0;
            arg2[2] = arg1Norm[0] * inverseLength;
            arg3[0] = arg1Norm[1] * arg2[2];
            arg3[1] = arg1Norm[2] * arg2[0] - arg1Norm[0] * arg2[2];
            arg3[2] = arg1Norm[1] * arg2[0];
        } else {
            double inverseLength = 1 / Math.sqrt(arg1Norm[1] * arg1Norm[1] + arg1Norm[2] * arg1Norm[2]);
            arg2[0] = 0;
            arg2[1] = arg1Norm[2] * inverseLength;
            arg2[2] = - arg1Norm[1] * inverseLength;
            arg3[0] = arg1Norm[1] * arg2[2] - arg1Norm[2]*arg2[1];
            arg3[1] = - arg1Norm[0] * arg2[2];
            arg3[2] = arg1Norm[0] * arg2[1];
        }
    }

    private static int computeRank(double[][] Matrix) {
        double abs, save, max = -1;
        int  maxRow = -1, maxCol = -1;
        for (int row = 0; row < SSA_LENGTH; row++) {
            for (int col = row; col < SSA_LENGTH; col++) {
                abs = Math.abs(Matrix[row][col]);
                if (abs > max) {
                    max = abs;
                    maxRow = row;
                    maxCol = col;
                }
            }
        }
        if (max < EPSILON) {
            return 0;
        }
        if (maxRow != 0) {
            for (int col = 0; col < SSA_LENGTH; col++) {
                save = Matrix[0][col];
                Matrix[0][col] = Matrix[maxRow][col];
                Matrix[maxRow][col] = save;
            }
        }

        double inverse = Matrix[0][maxCol];
        Matrix[0][0] /= inverse;
        Matrix[0][1] /= inverse;
        Matrix[0][2] /= inverse;

        if (maxCol == 0) {
            Matrix[1][1] -= Matrix[1][0]*Matrix[0][1];
            Matrix[1][2] -= Matrix[1][0]*Matrix[0][2];
            Matrix[2][1] -= Matrix[2][0]*Matrix[0][1];
            Matrix[2][2] -= Matrix[2][0]*Matrix[0][2];
            Matrix[1][0] = 0;
            Matrix[2][0] = 0;
        } else if (maxCol == 1) {
            Matrix[1][0] -= Matrix[1][1]*Matrix[0][0];
            Matrix[1][2] -= Matrix[1][1]*Matrix[0][2];
            Matrix[2][0] -= Matrix[2][1]*Matrix[0][0];
            Matrix[2][2] -= Matrix[2][1]*Matrix[0][2];
            Matrix[1][1] = 0;
            Matrix[2][1] = 0;
        } else {
            Matrix[1][0] -= Matrix[1][2]*Matrix[0][0];
            Matrix[1][1] -= Matrix[1][2]*Matrix[0][1];
            Matrix[2][0] -= Matrix[2][2]*Matrix[0][0];
            Matrix[2][1] -= Matrix[2][2]*Matrix[0][1];
            Matrix[1][2] = 0;
            Matrix[2][2] = 0;
        }

        max = -1;
        maxRow = -1;
        maxCol = -1;

        for (int row = 1; row < SSA_LENGTH; row++) {
            for (int col = 0; col < SSA_LENGTH; col++) {
                abs = Math.abs(Matrix[row][col]);
                if (abs > max) {
                    max = abs;
                    maxRow = row;
                    maxCol = col;
                }
            }
        }
        if (max <EPSILON) {
            return 1;
       }

       if(maxRow == 2) {
           for (int col = 0; col < SSA_LENGTH; col++) {
               save = Matrix[1][col];
               Matrix[1][col] = Matrix[2][col];
               Matrix[2][col] = save;
           }
       }

       inverse = 1 / Matrix[1][maxCol];
       Matrix[1][0] *= inverse;
       Matrix[1][1] *= inverse;
       Matrix[1][2] *= inverse;
       return 2;
    }

    private static double[][] computeEigenVectors(double[][] Matrix) {
        double[] eigenValues = computeEigenValues(Matrix);
        double[][] FirstValueMatrix = new double[SSA_LENGTH][SSA_LENGTH];
        double[][] SecondValueMatrix = new double[SSA_LENGTH][SSA_LENGTH];
        double[][] ThirdValueMatrix = new double[SSA_LENGTH][SSA_LENGTH];
        int rankFirst = 0;
        int rankSecond = 0;
        int rankThird = 0;
        double[][] result = new double[SSA_LENGTH][SSA_LENGTH];

        for (int row = 0; row < SSA_LENGTH; row++) {
            for (int col = 0; col < SSA_LENGTH; col++) {
                FirstValueMatrix[row][col] = Matrix[row][col];
                if (row == col) {
                    FirstValueMatrix[row][col] -= eigenValues[0];
                }
            }
        }

        rankFirst = computeRank(FirstValueMatrix);
        for (int row = 0; row < SSA_LENGTH; row++) {
            for (int col = 0; col < SSA_LENGTH; col++) {
                SecondValueMatrix[row][col] = Matrix[row][col];
                if (row == col) {
                        SecondValueMatrix[row][col] -= eigenValues[1];
                }
            }
        }

        for (int row = 0; row < SSA_LENGTH; row++) {
            for (int col = 0; col < SSA_LENGTH; col++) {
                ThirdValueMatrix[row][col] = Matrix[row][col];
                if (row == col) {
                    ThirdValueMatrix[row][col] -= eigenValues[2];
                }
            }
        }
        if (rankFirst == 0) {
            //M1 = A - evalue[1]*I;
           // M2 = A - evalue[2]*I;
            rankSecond = computeRank(SecondValueMatrix);
            rankThird = computeRank(ThirdValueMatrix);
            for(int i = 0; i < SSA_LENGTH; i++) {
                for(int j = 0; j < SSA_LENGTH; j++) {
                    result[i][j] = 0;
                    if (i == j) {
                        result[i][j] += 1;
                    }
                }
            }
            return result;
        }
        if(rankFirst == 1) {
            // evalue[0] = evalue[1] < evalue[2]
            rankSecond = computeRank(SecondValueMatrix);
            rankThird = computeRank(ThirdValueMatrix);
            double[] rowZero = Arrays.copyOf(FirstValueMatrix[0], SSA_LENGTH);
            double vectorNorm =  getVectorNorm(rowZero);
            for (double d: rowZero) {
                d = d / vectorNorm;
            }
            getSecondComplement(rowZero, result[0], result[1]);
            result[2] = computeCrossProduct(result[0], result[1]);
            return result;
        }
        //rankFirst == 2
        double[] rowZero = Arrays.copyOf(FirstValueMatrix[0], SSA_LENGTH);
        double vectorNorm =  getVectorNorm(rowZero);
        for (double d: rowZero) {
                d = d / vectorNorm;
        }
        //Obtainin Second Row - rowFirst
        double[] rowFirst = Arrays.copyOf(FirstValueMatrix[1], SSA_LENGTH);
        double vectorFirstNorm =  getVectorNorm(rowFirst);
            for (double d: rowFirst) {
                d = d / vectorFirstNorm;
        }
        //First Eigen Vector
        result[0] = getFirstComplement(rowZero, rowFirst);
        rankSecond = computeRank(SecondValueMatrix);
        if (rankSecond == 1) {
            rankThird = computeRank(ThirdValueMatrix);
            getSecondComplement(result[0], result[1], result[2]);
            return result;
        }
         //rankSecond == 2
        rowZero = Arrays.copyOf(FirstValueMatrix[1], SSA_LENGTH);
        vectorNorm =  getVectorNorm(rowZero);
        for (double d: rowZero) {
                d = d / vectorNorm;
            }
        //rowFirst
        rowFirst = Arrays.copyOf(FirstValueMatrix[1], SSA_LENGTH);
        vectorFirstNorm =  getVectorNorm(rowFirst);
            for (double d: rowFirst) {
                d = d / vectorFirstNorm;
            }
        //Second Eigen Vector;
        result[1] = getFirstComplement(rowZero, rowFirst);
        rankThird = computeRank(ThirdValueMatrix);
        result[2] = computeCrossProduct(result[0], result[1]);
        return result;
    }

    private double[][] computeRightEigenVectors( ) {
        double[][] svdMatrix = getSVDMatrix();
        double[][] hankelMatrix = getHankelMatrix();
        double eigenValues[] = computeEigenValues(svdMatrix);
        double leftEigenVectors[][] = computeEigenVectors(svdMatrix);
        final int TEMP_LENGTH = temperatureSeries.length;
        double[][] result = new double[SSA_LENGTH][TEMP_LENGTH - SSA_LENGTH + 1];
        for (int i = 0; i < SSA_LENGTH; i++) {
            for (int j = 0; j< TEMP_LENGTH - SSA_LENGTH + 1; j++) {
                result[i][j] = 0;
                for (int k = 0; k< SSA_LENGTH; k++) {
                    result[i][j] += hankelMatrix[i][k] * leftEigenVectors[j][k];
                }
                result[i][j] = result[i][j] / Math.sqrt(eigenValues[i]);
            }
        }
        return result;
    }

    public double[] extractAlphas() {
        double[][] svdMatrix = getSVDMatrix();
        double[][] eigenVectors = computeEigenVectors(svdMatrix);
        for (int i = 0;i < SSA_LENGTH; i++) {
            double norm = getVectorNorm(eigenVectors[i]);
            for (int j = 0; j < SSA_LENGTH; j++) {
                if (norm != 0) {
                    eigenVectors[i][j] = eigenVectors[i][j] / norm;
                }
            }
         }
        double lastCompNormSquare = 0;
        double[] result = new double[SSA_LENGTH-1];
        //v^2 = sum of last components
        for (int i = 0; i< SSA_LENGTH; i++) {
            lastCompNormSquare += eigenVectors[i][SSA_LENGTH - 1] *
                                    eigenVectors[i][SSA_LENGTH - 1];
        }
        //obtaining alphas
        for (int i = 0; i < SSA_LENGTH - 1; i++) {
            result[i] = 0;
            for (int j = 0; j < SSA_LENGTH; j++) {
                result[i] += eigenVectors[i][SSA_LENGTH - 1] * eigenVectors[i][j];
            }
            result[i] = result[i] / lastCompNormSquare;

        }
        return result;
    }
    private static double getTraceOfMatrix(double[][] arg) {
        double result = 0;
        for (int i = 0; i < arg.length; i++) {
            result += arg[i][i];
        }
        return result;
    }

    private static double getDeterminant(double[][]  arg) {
        double result = 0;
        result -= arg[0][2] * arg[1][1]* arg[2][0];
        result += arg[0][1] * arg[1][2] * arg[2][0];
        result += arg[0][2] * arg[1][0] * arg[2][1];
        result -= arg[0][0] * arg[1][2] * arg[2][1];
        result -= arg[0][1] * arg[1][0] * arg[2][2];
        result += arg[0][0] * arg[1][1] * arg[2][2];
        return result;
    }
    public static double[] computeEigenValues(double[][] arg) {
        double[] result = new double[SSA_LENGTH];
        double p1 = (arg[0][1] * arg[0][1]) + (arg[0][2] * arg[0][2]) + (arg[1][2] *arg[1][2]);
        //diagonal matrix
        if (p1 == 0) {
            result[0] = arg[0][0];
            result[1] = arg[1][1];
            result[2] = arg[2][2];
        } else {
            double q = getTraceOfMatrix(arg) / 3;
            double p2 = (arg[0][0] - q) * (arg[0][0] - q) + (arg[1][1] - q) * (arg[1][1] - q) +
                        (arg[2][2] - q) * (arg[2][2] - q) + 2 * p1;
            double p = Math.sqrt(p2 / 6);
            double[][] bMatrix = new double[SSA_LENGTH][SSA_LENGTH];
            for (int i = 0; i < SSA_LENGTH; i++) {
                for(int j = 0; j < SSA_LENGTH; j++) {
                    bMatrix[i][j] = arg[i][j];
                    if (i == j) {
                        bMatrix[i][j] -= q;
                    }
                    bMatrix[i][j] = bMatrix[i][j] / p;
                }
            }
            double r = getDeterminant(bMatrix) / 2;
            double phi = 0;
            if (r <= -1) {
                phi = Math.PI / 3;
            } else if (r >= 1) {
                phi = 0;
            } else {
                phi = Math.acos(r) / 3;
            }

            result[0] = q + 2 * p * Math.cos(phi);
            result[2] = q + 2 * p * Math.cos(phi + 2 * Math.PI / 3);
            result[1] = 3 * q - result[0] - result[2];
        }
        Arrays.sort(result);
        return result;
    }

    public double[] forecatsTemp() {
        if(temperatureSeries.length < MIN_FORECAST_SIZE) {
            throw new NotEnoughValuesToForecastException("Not Enough data for forecaset " + temperatureSeries.length);
        }
        double alphas[] =  extractAlphas();
        int  lastElement = temperatureSeries.length - 1;
        double[] forecast = new double[SSA_LENGTH];
        forecast[0] = alphas[0] * temperatureSeries[lastElement]
                        + alphas[1] * temperatureSeries[lastElement - 1];
        forecast[1] = alphas[0] * forecast[0] + alphas[1] * temperatureSeries[lastElement];
        forecast[2] = alphas[0] * forecast[1] + alphas[1] * forecast[0];
        return forecast;
    }
}
